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ABSTRACT 

We study charmonium using the standard relativistic formalism in the quenched ap- 
proximation, on a set of lattices with isotropic lattice spacings ranging from 0.1 to 
0.04 fm. We concentrate on the calculation of the hyperfine splitting between rjc and 
J/ijj, aiming for a controlled continuum extrapolation of this quantity. The splitting 
extracted from the non-perturbatively improved clover Dirac operator shows very 
little dependence on the lattice spacing for a < O.lfm. The dependence is much 
stronger for Wilson and tree-level improved clover operators, but they still yield 
consistent extrapolations if sufficiently fine lattices, a < 0.07 fm {aM{r]c) < 1), 
are used. Our result for the hyperfine splitting is 77(2) (6) MeV (where Sommer's 
parameter, ro, is used to fix the scale). This value remains about 30% below ex- 
periment. Dynamical fermions and OZI-forbidden diagrams both contribute to the 
remainder. Results for the rjc and J/ip wave functions are also presented. 



1 



1 Introduction 



Heavy QCD quarkonium systems have been thoroughly studied analytically within 
the heavy quark non-relativistic approximation, NRQCD and related heavy 
quark effective theories as pNRQCD |2j and vNRQCD [3 - for a recent review see jl]. 
These approaches, however, generically fail for charmonium. For cc the expansion 
parameter v^l(? of the effective theory is about 0.3, and higher order corrections in 
vjc seemingly become very large, even overwhelming the leading order terms. It is 
particularly challenging to reproduce the hyperfine splitting between the ^5*1 and the 
""^^o states, which for charmonium is M{J/il) — rjc) = 117 MeV. The lattice version 
of NRQCD predicts a value of the hyperfine splitting M{J/ip — rjc) = 55(5) MeV 
jH] far below the experimental value, although still quite remarkable taking into 
account all the approximations involved. Indeed lattice NRQCD, being an effective 
action with cutoff given by the heavy quark mass m^, does not allow for a continuum 
extrapolation: for the effective action to be valid one has to preserve arriq > 1. Given 
this, the estimation of the systematic uncertainties inherent in the discretization is 
quite difficult. This fact, together with the observation that the v/c expansion is 
not well justified for charmonium, encourages the suspicion that the origin of the 
discrepancy between experimental and computed hyperfine splitting could lie in the 
non-relativistic approximation. However, all other lattice determinations based on 
relativistic actions also underestimate the value of the hyperfine splitting, by as 
much as 30 — 50% jUl 13 El El HO] ■ This raises what one might call the puzzle of the 
hyperfine splitting. 

Almost all lattice calculations up to now, including NRQCD, have been per- 
formed within the quenched (valence quark) approximation (for recent reviews on 
heavy quarks on the lattice see [II1II21)- The effect of quenching has been estimated 
by looking at the predicted form of the hyperfine splitting in the heavy quark 
approximation 

M(W-,.) = 5?^|*(0)P , (1) 

with \l/(0) the value of the non-relativistic wave function at the origin. In Ref. JH] it 
is argued that the change in \I^(0) and in the running of due to dynamical quarks 
can give altogether a deviation of the quenched result from the real world case 
as large as 40%. This would make the hyperfine splitting a quantity particularly 
suitable for unveiling unquenching effects - typically the effects of quenching on 
other spectral quantities amount to only ~ 10%. However, first numerical results 
including dynamical quarks seem to indicate a much milder Nf dependence of the 
hyperfine splitting than what is needed to match the experimental result ^5 ■ 
This raises some worry about the reliability of the calculations performed up to 
now, in particular since they have been typically done at not so small values of 
the lattice spacing. This worry becomes more severe after noting that continuum 
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extrapolations of the hyperfine splitting seem to depend quite strongly on the kind of 
Dirac operator used for the calculation El El • Needless to say that the continuum 
limit is unique. Such strong dependence on the choice of Dirac operator has to reflect 
the existence of large lattice artifacts. 

Indeed, the reason why a reliable determination of the hyperfine splitting by 
means of purely non-perturbative relativistic calculations remains elusive is that 
charm is too heavy for most current lattice simulations. For charm the dominant 
lattice artifacts are of (9((amc)"), n > 1, with anic ~ 0.5 at the typical values 
of the lattice spacing. Lattice artifacts remain large and may completely spoil the 
determination of the charmonium spectrum. 

One particular approach that has been advocated in order to avoid large 0{amq) 
lattice artifacts for heavy quark systems is the use of anisotropic lattices [71. Such 
lattices ^Hj have different lattice spacings and at in spatial and temporal direc- 
tions, with at « as- In Ref. |7| it has been argued that by tuning appropriately the 
parameters of the action one could achieve reduced scaling violations of ©((aimq)"") 
while keeping still agrriq large. This is certainly a very attractive possibility, which 
would make anisotropic lattices quite advantageous over standard isotropic ones. 
In particular this is the approach that has been taken in several recent studies of 
the charmonium spectrum [7llHl 011011111 • However, recent analysis seem to point 
out that scaling violations governed by the large asTUq artifacts could revive at the 
one loop level fr7\ I18j. This casts doubt on the generic effectiveness of anisotropic 
actions in reducing lattice artifacts better than standard isotropic lattices. To set- 
tle this point, an analysis of the range of validity of the anisotropic approach, as 
performed in [2^, seems essential. 

Our methodology here is to compute M{J /ip — rjc) on very fine isotropic lattices 
within the quenched relativistic formalism. If the lattices used are sufficiently fine, 
this approach should allow for a controlled quenched continuum extrapolation, free 
of the systematic uncertainties which may have affected other determinations up to 
now. Particularly important for this will be the use of a non-perturbatively improved 
version of the Dirac operator, for which lattice artifacts are reduced to (9((amq)^). 
Preliminary results of this calculation have been presented in |21j . 

We want to stress that our calculation, which does not differ from any other 
in this respect, involves one other approximation besides quenching: OZI, meaning 
that Zweig-rule forbidden diagrams are not taken into account^. Such diagrams, 
as in Fig. Q contribute to the correlator of flavor singlet mesons like charmonium. 
They have been particularly studied for light quarks in connection with the U(1)a 
symmetry breaking and the 1]' mass, where the anomaly provides an enhancement of 
OZI contributions to pseudoscalar mesons - for lattice estimates see [221 123 • There 
is however no lattice calculation of their magnitude for heavy quarkonium states 

"'^We thank Stefan Sint for this remark. 
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like charmonium. In the quark model fI^ OH] these OZI amplitudes are expected 
to be proportional to the value of the wave function at the origin, and hence are 

— * 

suppressed for P-wave channels which have vanishing wave function at r = 0. For 
S-wave, their contribution is 0{al{mps)) and 0{al{mv)) suppressed for the pseudo- 
scalar and vector channels respectively. Experimental evidence of the goodness of 
the OZI rule and of this relative suppression comes from the hadronic widths of 
heavy pseudo-scalar and vector mesons. The effect of these diagrams is therefore 
expected to be small for heavy quark systems like charmonium. Here, we neglect 
such "disconnected" diagrams because of their high computer cost. We want to 
stress, however, that their contribution may amount to several MeV, perhaps a not 
so negligible effect given the small value of the hyperfine splitting. Moreover, the 
OZI contribution could be enhanced by mixing with glueballs with the appropriate 
quantum numbers provided these glueballs were almost degenerate in mass with 
the charmonium states, an effect that has been measured for light singlet mesons 

122 123- 




Figure 1: Zweig rule forbidden diagrams contributing to the pseudoscalar (Left) 
and vector (Right) channels. These closed quark loops are connected by gluonic 
contributions, of which we only show those that contribute to lowest order in ct^. 

The paper is organized as follows. Section 2 compiles our results. In Section 
2.1 details about our numerical simulations are presented, including the kind of 
Dirac operators used to measure meson masses and the strategy to extract the 
hyperfine splitting from 2-point correlators. In Section 2.2 we give our values for the 
^5*1— ^5*0 hyperfine splitting as well as the IP to'^S'i splittings, obtained with the non- 
perturbatively improved clover Dirac operator. Our continuum extrapolation of the 
charmonium spectrum is based on data at lattice spacings ranging from a = 0.0931 
fm to 0.0397 fm. A comparison with the results obtained from Wilson and tree-level 
clover Dirac operators is given in Section 2.3. In Section 2.4 we present the results 
for the rjc and J/ip wave functions. We investigate the volume dependence or our 
results in Section 2.5, by looking at the dependence on physical volume of both 
the spectrum and the wave functions. Finally, our summary and conclusions are 
presented in Section 3. 
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2 Results 



We have aimed at analyzing the S- and P-wave states indicated in Table ^ Pre- 
cise results will only be presented for the '^S'l — hyperfine sphtting, although 
preliminary results for P-wave meson masses will also be given. 



Table 1: Charmonium S- and P-wave states analyzed in this work. We indicate 
their quantum numbers and physical masses. Unless indicated experimental errors 
are one or less in the last indicated digit. 



Name 


(2s+l)^^ 


jPC 


Mass(GeV) 


Vc 




0-+ 


2.980(2) 


J/ij 




1 — 


3.097 


he 


'Pi 


1+- 


3.526 


XcO 


'Po 


0++ 


3.415 


Xcl 


'Pi 


1++ 


3.511 



2.1 Simulation details 

Gauge configurations are generated with the standard Wilson action. Simulation 
parameters - see Table |21 - have been chosen so as to fix a physical lattice size of 
about 1.6 fm. In addition we have also simulated a very fine 32'^ x 96 lattice at 
(3 = 6.6, corresponding to a smaller size of 1.3 fm (possible finite volume effects will 
be discussed in section The physical scale is set by tq from |77] . 



Table 2: Simulation parameters. The scale is set by tq from ^7\. The non- 
perturbative value of the clover coefficient, c^^ , has been computed in pHj . 





L^xT 


a(fm) 


La (fm) 


„NP 
^sw 


# conf. 


6.0 


18^ X 48 


0.0931 


1.68 


1.769 


190 


6.2 


24^ X 72 


0.0677 


1.62 


1.614 


90 


6.4 


32^ X 96 


0.0513 


1.64 


1.526 


60 


6.6 


32^ X 96 


0.0397 


1.27 


1.467 


130 



Quark propagators are computed using Wilson, tree-level clover and non- pertur- 
batively improved clover Dirac operators. The non-perturbative value of the clover 
coefficient is taken from |28j . 

A complete analysis of the spectrum, at all (3 values, has only been performed for 
the non-perturbatively improved clover Dirac operator. For Wilson and tree-level 
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clover Dirac operators we have results at /? = 6.2 (for Wilson), (3 = 6.4 and (3 = 6.6. 
The continuum extrapolations in such cases are based on available data at (3 = 6.0 
and 6.2 from UKQCD 0121. 
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Figure 2: Left: ^5*1 — ^Sq hyperfine splitting effective mass. Right: ^Pi — ^Si 
splitting effective mass. Data are obtained using the non-perturbatively improved 
clover Dirac operator. The dashed horizontal line is the experimental value and the 
band comprises our fit to a constant with its errors. 

We compute correlation functions of hadronic operators with zero-momentum 
point source and point sink. A time dependent effective mass M{t) is extracted by 
fitting the ratio of correlators at a pair of points {t, t + 1) to: 



(C(t)) _ cosh[(f - t) M{t)] 
(C(t + 1)) ~ cosh[(f M{t)] 



(2) 



Here (.) means summation over Yang-Mills configurations. The plateau in the ef- 
fective mass is then fitted to a constant which provides our estimate of the meson 
mass. 

In addition, we extract mass differences between the states in Table ^ from the 
ratio of correlators, which behaves as: 



RA{t) _ cosh[(f - t)MA{t)] cosh[(f - t - l)MB{t)] 
Rsit) ~ cosh[(f - t)MB{t)] cosh[(f - t - l)MA{t)] 



-AAf(i) 



(3) 
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if the correlators are dominated by a single pole. This allows to extract a time 
dependent effective mass for the splitting, AM(t) = Msit) — M^it), which is fitted 
to a constant in the range where the effective mass plateau sets in. The quality of 
the data for the hyperfine splitting is shown in Fig. |21 left. A similar plot for the 
^Pi — ^Si splitting is presented in Fig. [fright. Although the signal extracted from 
Eq. is less noisy than the one from Eq. (0), point sources are still too noisy for 
P-wave states and, within our statistics, do not allow a precise determination of the 
masses. 

2.2 Charmonium spectrum from non-perturbatively 
improved clover Dirac operator 

We present first the results obtained with the non-perturbatively improved clover 
Dirac operator, for which scaling violations in spectral quantities are expected to 
be 0{a^). A detailed comparison with Wilson and tree-level clover Dirac operators 
will be presented in section 12.31 

Table 3: Charmonium spectrum from non-perturbatively improved clover Dirac 
operator. Results are given in MeV with the scale set by tq from |27j. k has been 
tuned to maintain an approximately constant mass M{^Si) ~ 3095 MeV Va. The 
last two columns show our continuum extrapolation and the experimental value. 





(3 = 6.0 
K = 0.11865 


(3 = 6.2 
K = 0.12457 


(3 = 6.4 
K = 0.12755 


(3 = 6.6 
K = 0.12943 


a ^ 


Exp. 




3023(2) 


3019(3) 


3034(3) 


3014(3) 




2980 




3091(2) 


3093(3) 


3109(4) 


3085(3) 




3097 


3 c 1 C 
Oi — Do 


68.5(1.3) 


74.2(1.5) 


75.2(2.5) 


73.8 (2.1) 


77.2(1.7) 


117 


'Pi - 's, 


417(25) 


460(34) 


433(37) 


417(40) 


441(25) 


429 


'Po - 


342(19) 


352(25) 


369(34) 


397(41) 


387(14) 


318 


'Pi - 'S^ 


390(19) 


413(34) 


451(41) 


417(29) 


437(16) 


414 



We have used two different ways of setting the charm quark point: imposing 
that either roM(^S'o) or roM{^Si) equals the physical value. These two different 
choices allow to study the ambiguities in the scale determination induced by the 
approximations we have made: quenched and OZI. We want to stress again that our 
calculation is affected not only by quenching ambiguities but also by OZI effects, i.e. 
only a subset of the relevant quenched diagrams has been included. Taking this into 
account we expect the choice of J/ ip mass as reference to be a better one, since in that 
case OZI contributions to the J/ip mass are expected to be considerably suppressed 
compared to the pseudoscalar ones - see for instance and the discussion in 
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Section n For these reason the numbers presented here refer to the choice of J /if) 
mass as reference. The hyperfine sphtting with 71^(^5*0) fixed to the experimental 
Tjc value will be used below when we estimate the systematic uncertainty in our 
calculation. An alternative possibility, which we have not explored, would be to fix 
the charm quark point by using the Ds meson mass, as done for instance in heavy- 
light spectroscopy and for the determination of the charm quark mass [301 IHT] . This 
would be free of OZI ambiguities. 

Our results are collected in Table El As indicated above, the charm quark point 
has been obtained by fixing Mi^Si) ~ 3095 MeV for all a. We compute the hyperfine 
splitting from: {A) the difference between the '^Si and ^Sq masses and {B) the 
hyperfine splitting effective mass extracted from the ratio of correlators Eq. (jHj). 
Both determinations are perfectly consistent within errors but (B) is more precise 
and will be used henceforth. Fig. El shows the results for the '^5'i — ^Sq hyperfine 
splitting as a function of a^. The lattice spacing dependence is very small. On 
the coarsest lattice, for which aM{}So) ~ 1.4 and large scaling violations could be 
expected, the deviation from our continuum extrapolation amounts to only 11%. 
The cutoff dependence is well fitted linearly in a^. The continuum extrapolation 
is M{^Si - ^Sq) = 77(2) MeV. Excluding from the fit the point at /? = 6.0 gives 
M{^Si — ^Sq) = 74(1) MeV. If, instead of the vector mass, we fix the charm scale 
by setting M(^So) ~ 2.945 MeV Va, the result for the hyperfine sphtting goes up by 
6%. Including both these results as systematic error in our determination we quote 
as value of the hyperfine splitting from the non-perturbatively improved clover Dirac 
operator M(^Si — ^Sq) = 77(2)(6). This is about 30% below the experimental value. 

When comparing our number to previous lattice determinations - see 130 El 
- it is important to note that we have used tq to set the scale. It is is quite com- 
mon for charmonium analysis to use instead the spin averaged splitting, ^Pi — S 
or P — S, which is claimed to lead to considerably larger values of the hyperfine 
splitting. Variations with the scale input are usually blamed on the quenched ap- 
proximation ^. Although this is partially true, the large differences often quoted 
are also coming from large scaling violations. This has been already illustrated 
by the latest CP-PACS result P, obtained using the Fermilab anisotropic action. 
The discrepancy in the determination of the hyperfine splitting from ro or P — S 
has been reduced from 27 to 16% as the lattices used in the continuum extrapo- 
lation have changed from a e [0.099,0.208] fm to a e [0.0697,0.1374] fm 0. 
Their final number is M{^Si - ^So)p_s = 85.3(4.4)j+25j MeV. From ro they ob- 
tain instead M(^S'i — ^So)ro = 72.6(9)|^3'gj MeV (the previous result in |H1 was 
M{^Si — ^So)ro = 65(1) MeV, the discrepancy clearly reflecting the large systematic 
ambiguities due to scaling violations). Our final number, M(^Si — ^5'o)ro = 77(2) (6), 
lies in between their two determinations, and is compatible within errors with both. 

^Note, however, that this quantity is also afleeted by the OZI approximation 
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Figure 3: Continuum extrapolation of the hyperfine splitting with the non- 
perturbatively improved clover Dirac operator. The bare quark mass is tuned to 
maintain an approximately constant mass M(^S'i) ~ 3095 MeV Va. 
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Figure 4: Comparison between the experimental charmonium spectrum and our 
numerical data. The scale is set by tq. The value of the J/ip mass is fixed to the 
experimental one. 
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We have also computed the sphtting between P-wave states and the ^Si. Results 
are presented in Table |21 With point sources and within our limited statistics, our 
results are too noisy to attempt a reasonable continuum extrapolation and make a 
definite statement about scaling violations in these quantities. Still, our estimates 
of the J = 1 to J/ip splittings are consistent, within errors, with experiment. The 
discrepancy turns out to be larger for the ^Pq state. 

An estimate of the ambiguities inherent in the quenched plus OZI approximations 
can be extracted from the values of the continuum extrapolated spectrum. Figure 
E] presents a comparison between the experimental charmonium spectrum and our 
data. Deviations from experiment amount to 30% for the J/ip to t]c splitting, and 
22% for the ^Pq to J/ip splitting. They are considerably reduced for the J = 1, 
P-wave to J/ip splittings, which are consistent within errors with experiment. 

2.3 Dependence on the choice of Dirac operator and con- 
tinuum extrapolations. 

As discussed in the Introduction, it has been often reported El that the 
continuum extrapolated hyperfine splitting strongly depends on the choice of Dirac 
operator. This happens in particular with anisotropic actions like the ones used in 
[ZIIHIE] which are claimed to reduce OirUqasY scaling violations to O^rrigat)"', at -C 
tts- Since the continuum limit is unique, such strong dependence on the choice of 
Dirac operator necessarily refiects the existence of large lattice artifacts. Indeed, 
CP-PACS, in Ref. jH], concludes that m^ag < 1 is still certainly needed in order 
to obtain reliable continuum extrapolations with anisotropic lattices. This is not 
surprising, but it is an indication that anisotropic lattices might not, in general, 
succeed in removing scaling violations better than isotropic lattices do (indications 
that 0{g^mqas) indeed revive at the one-loop level have been reported in [T7| ITH]). 
As we will see next, even with the non-improved Wilson Dirac operator a reasonable 
estimate of the hyperfine splitting can be obtained if, but only if, lattices with 
spacing a < 0.07 fm, i.e. aMi^So) < 1, are used. 

The pseudoscalar mass and the hyperfine splitting from Wilson and tree-level 
clover Dirac operators are presented in Table El We only have data at /3 = 6.2 
(for Wilson), (3 = 6.4 and (3 = 6.6. To perform the continuum extrapolation we 
use /3 = 6.0 and 6.2 data from UKQCD Our strategy here has been slightly 

different than in the previous section. Instead of fine tuning the vector mass at each 
f3 value, we interpolate from a range of masses around the physical J/ip mass. We 
fit the dependence of M{^Si — ^Sq) versus l/M{^Si) at fixed P, and extract the 
value of the hyperfine splitting at the desired vector meson mass from the fit. This 
allows to extract the splitting at M(^5'i) = 3.095 GeV, the same vector meson mass 
which we fixed for the non-perturbatively improved Dirac operator. 

Concerning the continuum extrapolations we expect scaling violations to be: 
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Table 4: Pseudoscalar, vector mass and hyperfine splitting from Wilson and tree- 
level clover Dirac operator. Results are given in physical units (MeV) with the scale 
set by Tq. 





K /3 


Of) 




3e _ 1 Q„ 


\A/ 1 IcoTl 
V V llowli 




fi 9 


9728 ('3'! 


2764l'4'l 


35 ?(] 9) 


AA/i Icon 

V V lloVJll 




(\ 9 


2788 ('3') 


2821 (?,) 


33 2(0 ()) 


\A/ 1 IcoTl 
V V lloWii 




\j . ^ 


2Q1 3('3'l 


2943l'3'l 


30 4^0 


V V llOWll 


n 1350 


fi 2 


3oqqf4') 

Kjyj u u \ ^ I 


3125^4^ 


26 7f0 ()) 


V V llOVJll 


n 1 1RQ 


fi zL 


2925 f4') 


2963f4^ 


38 4 fl 6") 


\A/ 1 1 QOTl 

VV llovjli 


0.1380 


6.4 


3087('4'l 


31 99l'5'l 


35 4 n 7"! 


Wilson 


0.1371 


6.4 


3263f5) 


3288f4) 


28.5 (1.1) 


Wilson 


0.1415 


6.6 


2594(5) 


2652(6) 


60.4 (2.4) 


Wilson 


0.1400 


6.6 


2970(6) 


3017(7) 


48.3 (2.3) 


Wilson 


0.1385 


6.6 


3338(5) 


3375(5) 


39.1 (1.6) 


Wilson 


0.1375 


6.6 


3583(5) 


3613(5) 


32.3 (1.6) 


Clover 


0.1324 


6.4 


2931(4) 


2996(5) 


64.5 (2.8) 


Clover 


0.1320 


6.4 


3022(4) 


3081(4) 


59.2 (2.5) 


Clover 


0.1315 


6.4 


3129(4) 


3189(5) 


59.0 (2.4) 


Clover 


0.1335 


6.6 


2964(6) 


3031(6) 


68.9 (3.2) 


Clover 


0.1330 


6.6 


3116(6) 


3180(6) 


64.8 (3.0) 


Clover 


0.13225 


6.6 


3342(6) 


3399(6) 


58.9 (2.8) 



• 0{a) for the Wilson Dirac operator. 

• 0(0^) and 0{g'^a) for the tree-level clover Dirac operator. 

These are the dominant contributions but, for very coarse lattices, sub- leading terms 
might also be important especially since they depend on (am^)^. To test the ap- 
proach to the continuum, we present in Fig. El the results of two different extrapo- 
lations: 

(I) Top figure El using the three largest P values, 

• linearly in a for the Wilson Dirac operator 

• linearly in for the tree-level clover Dirac operator 

(II) Bottom figure El using all four j3 values, 
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Figure 5: Comparison between continuum extrapolation of M{'^Si — -^Sq) with 
Wilson, tree-level clover and non-perturbatively improved clover Dirac operators. 
Top: Wilson and tree-level clover improved data are fitted linearly in a and 
respectively. Bottom: Both Wilson and tree-level clover data are fitted with a plus 
dependence. Wilson data at = 6.0 and tree level clover data at (3 = 6.0 and 
p = 6.2 are based on UKQCD resuhs in 0121]. 
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• including a and terms for the Wilson Dirac operator 

• including a and terms for the tree-level clover Dirac operator 

We obtain as continuum extrapolations for the Wilson and tree-level clover re- 
spectively: Mi^Si-^So) = 67(6) MeV and 73(1) MeV from fit (I) and M{^Si-^So) = 
80(9) MeV and 99(5) MeV from fit (II). This is to be compared with the result from 
non-perturbative improvement M{^Si — ^Sq) = 77(2) MeV. 

Several remarks are in order here: 

(i) Only the non-perturbatively improved data show a weak dependence on the 
lattice spacing. Wilson and tree-level clover data are both significantly below their 
continuum extrapolations at all the simulated (3 values. 

[a) Linear extrapolations in a and for Wilson and tree-level clover respectively, 
including /3 = 6.0 (aM(^5'o) = 1.4) are not justified. 

{Hi) If too coarse lattices are included in the fit, continuum extrapolations be- 
come quite sensitive to the assumed dependence on the lattice spacing. This is 
clearly observed in the case of the tree-level clover data, for which fit (II) gives a 
considerably higher continuum value. It is only when the extrapolations start from 
sufficiently fine lattices that they come out reasonably consistent. A lower bound for 
consistent extrapolations seems to be a < 0.07 fm for Wilson and tree-level clover 
data, i.e. aM(^So) < 1. Surprisingly, though, non-perturbative clover improvement 
seems to work well even on the coarser (3 = 6.0 lattice where aM(^S'o) = 1.4. 

To summarize, non-perturbative clover improvement seems crucial to remove 
strong scaling violations. It is possible to extract comparable results from Wilson or 
tree-level clover improved Dirac operators, but very fine lattices are needed in order 
to remove the ambiguities inherent in the extrapolation procedure. 

2.4 Wave functions 

The origin of the the hyperfine splitting can be easily understood within the naive 
non-relativistic approximation (see for instance IM])- This approximation amounts 
to solving the Schrodinger equation in a non-relativistic Coulombic potential and 
dealing with relativistic corrections in perturbation theory. To zeroth order ^5*1 and 
^5*0 states are degenerate. The degeneracy is removed to first order in perturbation 
theory by the spin-spin interaction, giving a value of the hyperfine splitting: 



with \E'nr(0) the value of the non-relativistic wave function at the origin (\E'nr('") = 
(Sttp^)""*^/^ exp{— r/(2p)}, p = (AasiTig/S)^^). Perturbative corrections to the wave 
function also depend on the spin; to lowest order in perturbation theory, the value 
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of the wave function at the origin increases for the pseudoscalar and decreases for 
the vector according to |M1 

,1 



1+5np + (--^)^^|^)v1/nr(0) 



^nr(O). 



(5) 
(6) 



'6 ' 9 

where z/ ^ 7.241 x 10~^ and ^^(/i^) the strong running couphng evaluated at scale 
yU^. Here (5np denotes the, spin-independent, non-perturbative correction to the wave 
function at the origin (estimated in [85| l86]). 





J/tp wave function 



J I I I 



0.2 0.4 0.6 0.8 1 1.2 1.4 

r (fm) 



0.2 0.4 0.6 0.8 1 1.2 1.4 

r (fm) 



Figure 6: Scaling analysis of pseudoscalar (left) and vector (right) matter wave 
functions. The vertical scale is logarithmic. 

As in Ref. p8j, we have extracted gauge invariant wave functions from lattice 
4-point functions: 

|vl/(f)|2 = {ccM^,m{4i^,){r)\cc) = (7) 

Ar^(Tr { S{6,0;x,h) Q S{0,2ti;x,hyr^r, 

X 

5(0, 2ti; x + f, ti) Qt 5(0, 0; f + f, hY^/.T^} ) 

S is the quark propagator and T = {75,71} for {^So,'^Si}. We have measured 
matter and charge wave functions by setting Q = {1,70} respectively (SHI- A/" is a 
normalization constant fixed by imposing normalization 1 to the wave function. 
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Table 5: Parameters of the fit to Eq. (jSJ of the rj^ and J/ip matter wave function. 
We give both the direct values of the wave function at the origin ^'(O) and the values 
derived from the fit, ^1/(0), in GeV^. The parameter p is given in fermi. 



/? 


\%M\' 








6.0 
6.2 
6.4 
6.6 


0.156(6) 0.153(6) 1.24 (2) 0.155(4) 
0.124(7) 0.117(3) 1.24 (2) 0.172(4) 
0.108(5) 0.100(2) 1.20 (2) 0.176(3) 
0.117(9) 0.112(2) 1.14 (1) 0.164(2) 


0.093(4) 0.089(3) 1.46(3) 0.210(5) 
0.063(3) 0.060(1) 1.59(2) 0.255(2) 
0.045(3) 0.0445(3) 1.59(1) 0.283(1) 
0.049(3) 0.0500(5) 1.57(2) 0.273(2) 



A comparison between pseudoscalar and vector matter wave functions using the 
non-perturbatively improved clover Dirac operator is presented in Fig. IHl We have 
binned the wave function in bins of size a/2. We include in the figure results for 
(3 = 6.2, 6.4 and (3 = 6.6 which show very good scaling within errors, up to possible, 
small, finite size effects for (3 = 6.6 which will be discussed in the next section. 
The observed pattern corroborates qualitatively the predictions of the heavy-quark 
model: the value of the wave function at the origin increases for the pseudoscalar, 
decreases for the vector. Table El presents the results of the fit of matter wave 
functions {Q equal 1 in Eq. ((7|)) to the following ansatz motivated by the heavy 
quark non-relativistic approximation and the form of variational wave functions in 
potential models: 

|vi/(r)r = |^(0)rexp{-Q'}. (8) 

Our fits have been performed for r < 0.8 fm. In all cases we obtain a /ndf < 0.5 
(<1.4 for f3 = 6.0). In infinite volume, wave function normalization fixes the value 
of the wave function at the origin to |\l/(0)p = p/ (47rp^r[^]), a relation well satisfied 
by our fits. Given the good scaling properties of our /3 = 6.6 and 6.4 wave functions 
one can extract an estimate of the continuum values of the pseudoscalar and vector 
wave functions at the origin. Fitting to a constant the results for these two values 
of (3 we obtain \^r,M\^ = 0.110(4) and |^j/^(0)p = 0.047(2). This quantity is 
of phenomenological interest since it enters in many of the estimates of the heavy 
quark approximation as well as in potential models for heavy quarks. One example, 
apart from the hyperfine splitting, is the leptonic decay width of the J/ ip which, in 
the non-relativistic approximation, can be expressed as: 

r(J/^^e+e-) = l^^|vl/(0)p (9) 

with Qc the charm quark charge in units of the proton charge and a the fine struc- 
ture constant. Inserting our value of the J /if) wave function at the origin we obtain 
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T{J/ip e^e~) = 5.8 ± 0.2 keV, probably in too good agreement with the ex- 
perimental value TexpiJ/'ip — ^ e+e~) = 5.26 ± 0.37 keV, taking into account that 
formula does not include radiative corrections (appart from those affecting the 
wave function) which, estimated to lowest order, multiply the right hand side of Eq. 
© by a factor 1 - 16 a,(4m2)/(37r) ~ 0.5 [31E7|. 

We can also make use of formulas (0) and (jH)) to extract an estimate of the 
magnitude of non-perturbative contributions to the wave function at the origin (5np). 
For this we use mdmc) = 1.301(34) GeV from PU] and extract from the equations 
both ^NP and as- Plugging our values for |\E'^^(0)p and |\l/j/^(0)p gives 5np = —0.12 
and as = 0.761. The large value of ag needed to match our results is a clear 
indication that spin dependent, and hence relativistic, non-perturbative effects are 
indeed rather strong for these charmonium states. If we would instead fix the strong 
coupling constant to ag ~ 0.5 as in we would obtain a strong spin dependence 
of 5np which would moreover turn out to be of 0{1) for the rjc wave function. 

2.5 Finite volume effects 

Table 6: Pseudoscalar mass and hyperfine splitting from non-perturbatively im- 
proved clover Dirac operator. The lattice spacing is fixed to 0.093 fm (/? = 6.0) and 
the number of lattice points L, hence the physical volume La, is varied as indicated 
in the table. Results, averaged over 100 configurations (190 for L = 8), are given in 
physical units (MeV) with the scale set by vq. 



L 


La (fm) 






3 c 1 Q 


8 


0.75 


2958(10) 


3019(12) 


61.4(4.4) 


10 


0.93 


2953(5) 


3023(6) 


70.6(2.5) 


12 


1.12 


2957(3) 


3032(5) 


75.4(2.7) 


14 


1.30 


2947(3) 


3020(4) 


72.6(1.9) 


16 


1.49 


2952(3) 


3025(4) 


74.9(2.1) 


18 


1.68 


2949(2) 


3021(3) 


72.5(1.5) 



In this section we investigate how our results depend on the physical lattice 
volume. In particular we are interested in what happens with our finest (3 = 6.6 
lattice which has a somewhat small physical size La = 1.3 fm. An indication of the 
magnitude of finite volume effects can already be obtained from the wave function 
plots and fits in the previous section. Compare in Fig. IHl the results for (3 = 
6.6 (La=1.3 fm) and (3 = 6.4 (La=1.6 fm). Within errors, the wave functions 
do not show any clear signal of finite size effects, except for a small deviation of 
the pseudoscalar wave function at the tail, which is, however, not very significant 
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within our errors. One could also argue from the fits in Table that the wave 
functions at the origin are slightly larger than expected for the La=1.3 fm lattice, 
but this may be just a reflection, through the wave function normalization, of the 
finite volume effects at the tail. It is, anyway, clear from these plots that finite 
size effects on the charmonium wave function are really small on the La=1.3 fm 
lattice. In consequence, we do not expect the hyperfine splitting on this lattice to 
be significantly affected. 
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Figure 7: Finite volume effects on /3 = 6.0 pseudoscalar (left) and vector (right) 
wave functions. 



A complete finite volume analysis, including the volume dependence of the hyper- 
fine splitting in the continuum limit, is beyond the scope of this paper. To address 
the question of the relevance of finite volume effects, we have decided to study in- 
stead a set of L= 8, 10, 12, 14, 16 and 18 lattices at the coarsest, (3 = 6.0, lattice 
spacing. Table IHl presents our results for the pseudoscalar and vector masses and the 
hyperfine splitting for k = 0.11925, corresponding, for L=18, to a pseudoscalar mass 
of M{^Sq) = 2949(2)MeV. Our results are not precise enough to attempt a fit of 
the volume dependence. Still, the variation of the hyperfine splitting amounts to at 
most 4% for La > 1 fm. Further indication about the smallness of finite volume ef- 
fects for lattices larger than 1 fm comes from the study of the volume dependence of 
wave functions. Fig. [7| shows the results for the set of lattices in Table Deviations 
from the large volume behavior are negligible for La > 1.12 fm. 
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3 Conclusions 



Table [7| compiles our result for the J/ip—rjc liyperfine splitting together with previous 
determinations using other lattice formalisms [HI E] • Our result is obtained with the 
non-perturbatively improved clover Dirac operator. The systematic error covers the 
difference between choosing as reference scale for charm the rjc or J/ ip mass as well 
as different continuum extrapolations (including or not (3 = 6.0). Our final result is 
AM(J/^ - r],) = 77(2)(6)MeV. 



Table 7: Hyperfine splitting in MeV from different lattice approaches. The scale is 
set by vq or by the spin-averaged P — S splitting. The result quoted as 'This work' 
is obtained from the non-perturbatively improved clover Dirac operator. 





This work 


CP-PACS [9^ 


CP-PACS iHj 


latest NRQCDjni 


Formalism 


Relativistic 


Relativistic 


Relativistic 


Non relativistic 


Lattice 


Isotropic 


Anisotropic 


Anisotropic 




Extrapolation 


Continuum 


Continuum 


Continuum 




Scale 


To 




P-S 


P-S 


AM 


77(2) (6) 


72.6(9)j!5:^{ 


85.3(4.4)[!^:j; 


55(5) 



Our value for the quenched hyperfine splitting within the OZI approximation 
remains 30% below the experimental result. Dynamical quark effects are usually 
expected to be (9(10 — 20)%, which amounts for a large part but not all of the 
discrepancy. Actual dynamical quark effects may turn out to be larger for this 
particular quantity but the remaining discrepancy might also be due to the OZI 
approximation. 

Our study shows the virtue of the "brute force" approach: reliable, consis- 
tent continuum extrapolations can be obtained from improved or non-improved 
discretizations if the lattice is fine enough. Conversely, using a coarse lattice im- 
plies continuum extrapolations where non-leading terms may be significant, thus 
introducing a systematic error which is very hard to control. From our study, the 
boundary between these two regimes sits where one would expect, around aMgg ~ 1. 

Remarkably however, the non-perturbatively improved clover Dirac operator ap- 
pears to give reliable extrapolations even starting from aMgg ~ 1.4. Non-leading 
terms remain small, perhaps because leading corrections 0{a'^) are themselves very 
small. We consider our Fig. |31 as a spectacular advertisement for using this Dirac 
discretization. 

We believe our result finally closes the long debate on the magnitude of the 
quenched charmonium hyperfine splitting (within the OZI approximation). While 
NRQCD or an anisotropic discretization may yield a similar value, the reliability 
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of such a result remains questionable: in the first case, one deals with an effective 
theory where the continuum limit cannot be taken; in the latter, the advantage over 
isotropic lattices in removing 0{amq)"' scaling violations remains unclear. Until this 
point is settled, and although very accurate results can be obtained with anisotropic 
actions, this reduction in statistical errors can be more than offset by an increase in 
systematic errors due to possible scaling violations. Anisotropic lattices may be very 
useful in other contexts of course like for instance at finite temperature. Perhaps 
our result can be used to fine-tune the anisotropic actions involved. 

Note finally that the charmonium system is ideally suited for a precision study 
on the current generation of PC clusters. Large, quenched lattices can be simulated 
efficiently with little inter-node communication, and high accuracy can be obtained 
at small cost. This computer environment seems well suited for a measurement of 
OZI-effects. 
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